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State of the art variable composition structure prediction based on density functional theory demonstrates 
that two new stoichiometries of PN, PN 3 and PN 2 , become viable at high pressure. PN 3 has a 
skutterudite-like Immm structure and is metastable with positive phonon frequencies at pressures between 
10 and 100 GPa. PN 3 is metallic and is the first reported nitrogen-based skutterudite. Its metallicity arises 
from nitrogen p-states which delocalise across N 4 rings characteristic of skutterudites, and it becomes a good 
electron-phonon superconductor at 10 GPa, with a T c of around 18 K. The superconductivity arises from 
strongly enhanced electron-phonon coupling at lower pressures, originating primarily from soft collective 
P-N phonon modes. The PN 2 phase is an insulator with P2/m symmetry and is stable at pressures in excess of 
200 GPa. 



High pressure synthesis offers a route to new materials with promising technological applications such as 
superconductors 1 " 4 as well as to compounds with unexpected stoichiometries 5 . New theoretical structure 
prediction techniques 6 " 9 have been used to identify novel superconductors by applying pressure to stoi- 
chiometries that form phases without the properties of superconductivity, metallicity or even stability at atmo- 
spheric pressure 10 " 15 . In spite of their great potential, nitride-based materials remain relatively unexplored due to 
the high stability of nitrogen molecules at ambient pressure; the strong N=N triple bond results in a high kinetic 
barrier to polymerisation, thus extreme temperatures and pressures are required to synthesise many nitrides. 
However, at high pressures, the physical and chemical properties of nitrogen change substantially, allowing the 
formation of new polymeric crystalline phases 16 and the possibility of forming new nitrogen-based superconduct- 
ing crystals. 

Nitride-based superconductors were among the first discovered, and include examples with remarkably high 
superconducting critical temperatures (T c ), such as NbN (14.7 K) 17 and /?-HfNCl (25.5 K) 18 . However, in both 
instances, the superconductivity does not originate from the nitrogen atoms 19 , and the XNC1 family of layered 
compounds have been proposed to be unconventional superconductors 20 . 

Another important example of superconductivity in pnictogen-based compounds involves the skutterudites, 
though, as far as we know, no nitrogen-based skutterudite has been found yet. They have the empirical formula 
TX 3 , where T is typically a transition metal and X is a pnictogen. They form a low density body-centred cubic Im3 
lattice consisting of dodecahedral cages (Fig. 1). This cage can be filled with a guest atom to form 
MT 4 X 12 , potentially changing the properties of the material, making it semiconducting 21 , paramagnetic 22 ,super- 
conducting 23 or even thermoelectric 24 .They are also characterised by square X 4 4 ~ groups connecting the cages, 
which are thought to influence electrical conductivity and superconductivity 25 . Filled skutterudite superconduct- 
ing materials have been synthesised under pressure 26 , and found to have high T c values, for example LaRu 4 Asi 2 
with T c = 10.3 K 26 and La 0 . 6 Rh 4 Pi 2 with T c = 17 K, the highest known among metal phosphides 27 . It has been 
shown that although the guest atom in filled skutterudites stabilises the lattice, the superconductivity is an 
intrinsic property of the T 4 X 12 framework 28,29 . 
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Figure 1 | The skutterudite-like Immm structure of PN 3 in its body- 
centred orthorhombic conventional cell. The skutterudite structure has 
Im3 symmetry, but the PN 3 phase is characterised by an orthorhombic 
distortion. Phosphorus atoms are six-coordinated in both cases. 
Phosphorus atoms are yellow and nitrogen atoms purple. 

In this article we perform a systematic search for thermodynami- 
cally stable phosphorus nitrides under pressure using evolutionary 
structure prediction and density functional theory (DFT). Beyond the 
only known stable phosphorus nitride, P 3 N 5 , which forms a variety of 
pressure- dependent phases 30 " 33 , we find that a novel phosphorus 
nitride skutterudite (PN 3 ) containing six-coordinated phosphorus 
can be formed at undemanding pressures (—10 GPa). PN 3 is the first 
reported nitrogen-based skutterudite. It is metallic, metastable and 
dynamically stable between 10 and 100 GPa. Although it is unstable 
at ambient pressure, it may be possible to stabilise it using off- 
stoichiometry doping with a guest species at the centre of the cage, 
as with other skutterudite superconductors. PN 3 is a good supercon- 
ductor with a T c of around 18 K at 10 GPa, higher than for any 
known skutterudite. We also report the discovery of a new insulating 
phase with PN 2 stoichiometry, which is stable above 200 GPa. 
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Figure 2 | Enthalpy per formula unit as a function of pressure for the best 
phases with the P 3 N 5 stoichiometry. 

Results 

High pressure phosphorus nitride structures. In the first instance, 
fixed composition searches were performed using one and two 
formula units of P 3 N 5 at pressures of 25 GPa, 50 GPa, 100 GPa 
and 200 GPa, yielding the Imml y-phase at lower pressures, and 
the previously predicted C2/c oxyvanite phase at high pressures, in 
agreement with the findings of Dong et al 34 The recovery of an 
experimentally verified phase in addition to a structure predicted 
through chemical intuition is a reassuring vindication of the 
USPEX code. According to the results presented in Fig. 2, the a 
phase, which is stable at ambient pressure, undergoes a transition 
to the y phase at around 7 GPa, which transforms to the oxyvanite 
phase at 34 GPa. The previously predicted 3 ' phase 35 is at least 0.2 eV 
higher in enthalpy than y. 

Variable composition searches were performed at 100 GPa in 
order to identify potential stable stoichiometries; these were followed 
by fixed composition searches for the best stoichiometries at pres- 
sures of 25 GPa, 50 GPa, 100 GPa and 200 GPa using PN 3 , PN 2 , 
P 3 N 4 , P 3 N 5 and PN. The resulting convex hull is shown in the 




Figure 3 | Convex hull for the P x Nx_ x system at 100 GPa (top) and at 
different pressures (bottom). H is the enthalpy per atom of the generated 
structure, H[P] the enthalpy per atom of the stable phosphorus phase 
(simple hexagonal phosphorus at 100 GPa) and H[N] is the enthalpy per 
atom of the stable nitrogen phase (cubic gauche nitrogen at 100 GPa). The 
lower panel only includes the Jmmm-PN 3 , P2/ra-PN 2 and oxyvanite-P 3 N 5 
structures. 
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Figure 4 | The Pllm structure of PN 2 > which becomes stable at pressures 
of above 200 GPa. 

upper panel of Fig. 3, illustrating the most stable compositions: PN 3 , 
PN 2 and P 3 N 5 . The enthalpies of the most interesting structures are 
plotted at different pressures in the lower panel. A new insulating PI I 
m phase of PN 2 (structure in Fig. 4 and electronic density of states in 
Fig. 5) becomes stable with respect to N and P 3 N 5 at pressures in 
excess of 200 GPa, while the best PN 3 structure {Immm) is very close 
to the hull below 100 GPa. It is noteworthy that both new structures 
consist of octahedral PN 6 units, containing octahedral six-coordi- 
nated phosphorus atoms; this was unprecedented until a pressure 
induced transition to six-coordinated phosphorus was observed in 
A1P0 4 as a densification mechanism 36 . The structural parameters for 
these two new structures are listed in Table I. 

The Jrarara-PN 3 skutterudite. The skutterudite-like PN 3 is formed 
from dodecahedral cages (Fig. reffig: Immm- structure) consisting of 
six- coordinated phosphorus atoms at the vertices of a cuboid, with 
the faces occupied by N 2 pairs, forming pentagonal rings. In contrast 
with all other known and predicted PN structure, the nitrogen atoms 
are all three-coordinated, compared with the y, oxyvanite and new 
P2/m phases, which contain both three- and four-coordinated 
nitrogens. Surprisingly, the N-N bonds are extremely short; at 
100 GPa, they are 1.28 A, increasing to 1.30 A at 10 GPa. No other 
structures were found with such a short bond length. As the pressure 
decreases, the orthorhombic distortion in the Immm cell becomes 
more pronounced in the b direction, causing the structure to break 
down; it is likely that this is thermodynamically favourable since the 
N 2 units are very close in geometry to their configuration in the gas 
phase. However, the structure is dynamically stable with well-defined 
positive phonons at pressures as low as 10 GPa. Thus the possibility 
remains that the skutterudite structure may be synthesised at 
moderate pressures. 

MX 3 skutterudites are also characterised by six square X 4 rings in 
the conventional 32-atom cell. In PN 3 , these rings are rectangular 
rather than square as a result of the orthorhombic distortion, but can 
be discerned as pairs of N 2 units separated by 1.9 A at 100 GPa, 
which could be either a weak bond or a non-bonding contact. A 
Bader analysis reveals that the nitrogens have little ionic character, 
with a charge of 0.065 to 0.075e. These X 4 rings have been shown to 
give rise to conduction and superconductivity in other skutterudites, 
namely CoP 3 and NiP 3 37 and LaFe 4 P 12 25 . 

PN 3 is metallic throughout its stable pressure range. There is a 
distinct set of three bands (Fig. 6) dispersing from approximately 
1 eV below e F to 2 eV above e F responsible for the conductivity. The 
partial density of states (Fig. 6) demonstrates that the dominant 




Energy (eV) 

Figure 5 | Electronic density of states of the Pllm phase at 200 GPa. 

contribution to the conducting bands arises from nitrogen 2p atomic 
orbitals. It is interesting to note that in electronic structure calcula- 
tions of other metallic skutterudites, namely CoP 3 and NiP 3 , the 
highest occupied bands consist of phosphorus centred p-orbitals 37 , 
which occupy analogous positions to nitrogen in PN 3 , phosphorus 
being isoelectronic to nitrogen. 

The three bands around the Fermi level were projected onto max- 
imally localised Wannier functions (MLWF), as illustrated in Fig. 7. 
The Wannier projections have p- orbital -like symmetry, but are delo- 
calised across N 2 pairs which form adjacent skutterudite N 4 rings, 
accounting for the conductivity. 

Phonon spectra and electron-phonon coupling in the Iraram-PN 3 
skutterudite. The phonon spectra of PN 3 (Fig. 8) are characterised 
by two groups of modes; firstly dispersive low-energy modes that 
involve collective vibrations of P and N atoms. These modes are 
separated from the second group, high-energy Einstein-like modes, 
primarily corresponding to vibrations of the short N—N bonds, by a 
distinct gap. As the pressure is lowered, the acoustic modes and some 
low- energy optical modes related to distortions of the whole 
dodecahedral cage are softened. This gives rise to a structural 
instability at ambient pressure as the lowest energy mode at the X 
and R point acquire an imaginary frequency. The instability is related 
to the decomposition of PN 3 and the formation of N 2 pairs. 
Structurally, this is manifested as the formation of PN "sheets" 
enclosing N 2 molecules, and is reflected by significant changes in 
the phonon dispersion. However, at 10 GPa the phonon spectrum 
does not display any imaginary frequencies and the crystal structure 
is dynamically stable. 

The electron-phonon coupling constant k, calculated from equa- 
tion (3), scales linearly with the density of states at the Fermi level and 
is inversely proportional to the square of the phonon frequencies. Its 
calculated values are presented in Table II and the Eliashberg func- 
tions are plotted in Fig. 9. At all pressures studied, the eDOS at e F 
remains high due to the presence of weakly dispersive bands in the 
vicinity of the Fermi level (Fig. 6). This ensures a large electronic 
contribution to the Eliashberg functions coming from the double 
electronic Dirac delta in equation (1). Thus, the evolution of k under 
pressure is mostly determined by the evolution of the phonon 
spectra. 

At 10 GPa the electron-phonon coupling constant reaches a 
remarkable value of 1.05. The low-energy modes (<400 cm" 1 ) con- 
tribute more strongly to k due to their low frequencies. This is evi- 
denced by the strong contribution of the low- energy part of the 
Eliashberg function to k, demonstrated in Fig. 9. The contribution 
to k of the lowest energy acoustic mode at the R special point is 
particularly large. At 20 GPa, the electron-phonon coupling 
decreases as low-frequency modes gain energy. This causes the 
strong peaks in ot 2 F(w) at low energy to shift to higher frequencies 
and thus reduces their contribution to k. At higher pressures, the 
energy of the lowest energy acoustic mode at R is strongly increased 
and, consequently, its contribution to k is strongly suppressed, indu- 
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Table 1 

Phase 


Structural parameters 
Pressure (GPa) 


for the two new PN phases 
Lattice a, b, c (A) a, p, y (°) 


Wyckoff position 




Atomic coordinates (fractional) 
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0.0000 
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0.0000 
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0.1778 


0.5000 


0.6075 



cing a considerable reduction of the electron-phonon coupling 
constant. 

Using the computed values of the electron-phonon coupling con- 
stant, we estimated the superconducting T c using the standard values 
of /n* = 0.10 and [i* = 0.13 for the Coulomb pseudopotential in 
equation (2). These values are summarised in Table II. T c reaches 
18.6 K for /a* = 0.10 at 10 GPa, a high value for an electron-phonon 
coupling superconductor. Using fi* = 0.13, a slightly lower but sim- 
ilar 16.1 K value is obtained, showing that for large values of X T c is 
weakly dependent on [i*. The evolution of T c is consistent with the 
pressure dependence of X and both are suppressed at high pressures. 

Inclusion of a guest atom in the skutterudite lattice. Filled skut- 
terudites are of particular interest for their potential thermoelectric 
properties. If the void in a skutterudite cage is significantly larger 
than the guest atom, then the vibrational modes of the guest atom are 
expected to be strongly anharmonic, resulting in "rattling modes" 24 . 
These modes tend to reduce thermal conductivity by scattering the 
phonons responsible for heat transport 38 and, consequently, enhance 
the figure of merit of thermoelectric materials 39 . We now turn to the 
question of whether the skutterudite remains stable when it hosts 
such a guest atom. 




It should be noted that the cage in PN 3 is considerably smaller than 
for other skutterudites, which are large enough to hold atoms such as 
barium, which has a covalent radius of 2.53 A. PN 3 has a body- 
centred orthorhombic lattice at lower pressures, in contrast with 
most skutterudites which are body-centred cubic; the distance 
between the body centre of the cage and the closest point on the 
surface of the cage (specifically, the midpoint of a N 2 bond) is 1.81 
A, compared with 2.79 A for CoAs 3 . This greatly limits the choice of 
guest atom to smaller species such as hydrogen, lithium and 
beryllium. 

A beryllium atom proved to be too large, resulting in imaginary 
frequencies throughout the Brillouin zone. A lone hydrogen atom 
was too thermodynamically unstable to occupy the body centre of the 
lattice, drifting towards a phosphorus atom during geometry opti- 
misation and breaking the Immm symmetry. The phonon dispersion 
for PN 3 with a lithium guest has a soft mode at the X special point, 
and is therefore not stable; this was the case through the pressure 
range 10-100 GPa. In conclusion, it seems there is no stable XP 4 N 12 
filled skutterudite. 

The inclusion of a guest atom and a reduction of the external 
pressure have a similar effect on the atomic structure: the guest atom 
forces the cage to expand, just as the structure loses its dynamical 
stability as the lattice parameter increases at pressures below 10 GPa. 
It maybe possible to stabilise the structure by doping the system with, 
for example, one lithium atom per eight conventional cells (i.e. a 2 X 
2X2 supercell), although this would necessitate phonon calculations 
that would be prohibitively expensive. 

Discussion 

In conclusion, we used density functional theory coupled with evolu- 
tionary structure prediction to explore the variable composition con- 
figuration space of phosphorus nitride. We identified a new stable 
stoichiometry, insulating PN 2 with P2/m symmetry, at pressures in 
excess of 200 GPa, and a fascinating new metastable structure, PN 3 , 
which is dynamically stable between approximately 10 and 100 GPa. 

The PN 3 structure has Immm symmetry, belongs to the Im3 skut- 
terudite class of compounds, and is the first to contain nitrogen as 



Figure 6 | Electronic structure of the Immm PN 3 phase. The electronic 
density of states (top panel) is decomposed into contributions from P 
atoms, N 2s electons and N 2p electrons. The band structure near the Fermi 
level is displayed at 10 GPa, 20 GPa, 50 GPa and 100 GPa. 




Figure 7 | A MLWF projection of the states from 1 eV below to 2 eV 
above e F at 100 GPa The orbital is spread over skutterudite N 4 rings, 
giving rise to electronic conductivity. 
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Figure 8 | Phonon dispersion and phonon density of states for the Immm phase at 0 GPa, 10 GPa and 100 GPa. 



pnictogen. It is characterised by an orthorhombic distortion, the 
extent of which increases as the pressure is decreased. Skutterudite- 
PN 3 is metallic throughout the range of pressures examined, 10- 
100 GPa, its electrical conductivity arising from the delocalisation 
of nitrogen p-states across N 4 rings. Most significantly, it becomes a 
good electron-phonon superconductor with a T c of around 18 K at 
10 GPa due to strong electron-phonon coupling with low-frequency 
phonon modes. At higher pressures the superconducting T c is sup- 
pressed as low-energy acoustic modes gain energy. The supercon- 
ducting behavior of PN 3 is comparable to other high-pressure 
superconductors such as Te, Li and S, where the evolution of T c under 



Table II | Calculated k, a>| og and T c values for the /mmm-PN3 skut- 
terudite. T c is calculated with two different values for the Coulomb 
pseudopotential fi* 



P (GPa) 


X 


co| og (meV) 


T C [^ = 0.10](K) 


rj/i* = 0.13](K) 


10 


1.05 


21.5 


18.6 


16.1 


20 


0.77 


31.4 


15.5 


12.4 


50 


0.41 


56.4 


3.0 


1.4 


100 


0.39 


59.9 


2.3 


0.9 



pressure is determined by the behavior of low-energy phonon modes, 
and is enhanced close to dynamical instabilities 40 " 42 . 

The phonon dispersion at 10 GPa shows no sign of instability, 
opening the door to the possibility of its synthesis, perhaps by chem- 
ical synthesis from precursors in the correct ratio, or thin film depos- 
ition on a substrate with a similar lattice parameter. Furthermore, 
although the dodecahedral cage is too small to host one guest atom 
(H, Li, Be) per primitive cell, it may be possible to stabilise the PN 3 
lattice at ambient pressure using off-stoichiometry doping with a 
guest species. 

Methods 

The USPEX code 6,7 was used in conjunction with DFT to determine the most ther- 
modynamically stable PN crystal structures for several different compositions. 
USPEX employs evolutionary algorithms to find the atomic structure with the lowest 
enthalpy at a given pressure, given knowledge only of the chemical composition of the 
unit cell. It has been employed with great success with many other systems at high 
pressure 43 " 46 . 

All structures were optimised with the CASTEP DFT code 47 at the GGA level, using 
the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 48 . All lattice 
parameters and atomic positions were relaxed without constraints until forces varied 
by less than 0.05 eV A" 1 and energies by less than 1 X 10" 5 eV between BFGS steps. 
SCF steps were considered to be converged when the total energy varied by less than 1 
X 10" 6 eV between steps. Electronic wavefunctions were expanded using a plane 
wave basis set with a cutoff of 350 eV, and core electrons (Is 2 in the case of N and 
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Figure 9 | Eliashberg function a l F(co) and the integrated electron-phonon coupling constant X{co) 
10 GPa, 20 GPa, 50 GPa and 100 GPa. 



da/a 2 F(a/) / co' for the Immm phase at 



ls 2 2s 2 2p 6 in the case of P) were represented with ultrasoft pseudopotentials. 
Monkhorst-Pack grids were generated per cell, using a spacing of In X 0.05 A -1 . 
Ultrasoft pseudopotentials were employed with the intention of minimising cal- 
culation times; the nitrogen potential in particular had a core radius of 1.4 a 0 , which 
may have resulted in core overlap in the case of structures with extremely short N— N 
bonds; however, enthalpies were compared with a hard reference pseudopotential 
with a 1.0 a 0 core radius, and the results found to be consistent. 

Phonon and electron-phonon calculations were performed on the 16-atom 
primitive cell of the Immm skutterudite structure using density- functional per- 
turbation theory (DFPT) 49 as implemented in the Quantum Espresso code 50 . 
Ultrasoft pseudopotentials, a plane wave cutoff of 60 Ry, a density cutoff of 
600 Ry and a 4 X 4 X 4 Monkhorst-Pack grid for the electronic integrations 
were used throughout. The cell was fully reoptimised using the new pseudo- 
potentials, and selected conduction bands were projected onto maximally loca- 
lised Wannier functions (MLWF) using the wannier90 code 51 ' 52 . Phonon 
frequencies and electron-phonon coefficients were calculated on a 2 X 2 X 2 q- 
point grid. 

The Eliashberg functions of the electron-phonon coupling were calculated as 



a 2 F(w) -- 



N(0)N k N q 



1 Itqnmfi 

x5(co-a) (1 (q)). 



J2 k,,k +q J <5(ek„)<5(e k + qm) 



(1) 



In % (1) Sk„,k + q „, = E» M5V/5u"(q)\k + q m)^(q) / (q) is the 
electron-phonon matrix element, where \kn) is a Kohn-Sham state with energy e^ n 
measured from the Fermi level (e f ), N k and N q are the number of electron and phonon 
momentum points used for the Brillouin zone sampling, N(0) is the density of states 
per spin at e F , Vis the self-consistent Kohn-Sham potential, u s *(q) is the Fourier 
transformed displacement from equilibrium of atom s along Cartesian direction a, M s 
is the mass of atom s, and co^Cq) and e^ a (q) are respectively phonon frequency and 
polarisations with momentum q. A finer 12X 12X 12 mesh was used in the sum over 
k-points in Eq. (1). Superconducting T c values were estimated from the Allen-Dynes 
modified McMillan equation 53 : 



T c - 



a»i 0 g 



exp 



1.04(1+ A) 
k-li* (1-0.62/1) 



where the electron-phonon coupling constant is calculated as 

(x 2 F(co) 



1 = 2 dco- 



(2) 



(3) 



the logarithmic frequency average is a»i og = exp 0 J dco 
Coulomb pseudopotential parameter. 



a 2 F(co) 



In (D ) , and fj* is the 
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